Diversity and noise effects in a model of homeostatic regulation of the sleep- wake cycle 
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Recent advances in sleep neurobiology have allowed development of physiologically based mathe- 
matical models of sleep regulation that account for the neuronal dynamics responsible for the regu- 
lation of sleep-wake cycles and allow detailed examination of the underlying mechanisms. Neuronal 
systems in general, and those involved in sleep regulation in particular, are noisy and heterogeneous 
by their nature. It has been shown in various systems that certain levels of noise and diversity can 
significantly improve signal encoding. However, these phenomena, especially the effects of diversity, 
are rarely considered in the models of sleep regulation. The present paper is focused on a neuron- 
based physiologically motivated model of sleep-wake cycles that proposes a novel mechanism of the 
homeostatic regulation of sleep based on the dynamics of a wake-promoting neuropeptide orexin. 
Here this model is generalized by the introduction of intrinsic diversity and noise in the orexin- 
producing neurons, in order to study the effect of their presence on the sleep- wake cycle. A simple 
quantitative measure of the quality of a sleep-wake cycle is introduced and used to systematically 
study the generalized model for different levels of noise and diversity. The model is shown to exhibit 
a clear diversity-induced resonance: that is, the best wake-sleep cycle turns out to correspond to 
an intermediate level of diversity at the synapses of the orexin-producing neurons. On the other 
hand, only a mild evidence of stochastic resonance is found, when the level of noise is varied. These 
results show that disorder, especially in the form of quenched diversity, can be a key-element for an 
efficient or optimal functioning of the homeostatic regulation of the sleep-wake cycle. Furthermore, 
this study provides an example of a constructive role of diversity in a neuronal system that can be 
extended beyond the system studied here. 
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All biological systems are inherently noisy and hetero- 
geneous. Disorder is mostly expected to disturb proper 
functioning of a system, like it can be the case with noise 
in a radio signal. However, it has been demonstrated 
by numerous studies that noise can actually improve sig- 
nal encoding - the so-called stochastic resonance phe- 
nomenon. Recently, it was discovered that quenched di- 
versity (heterogeneity) can also enhance the response of a 
system to an external perturbation (diversity-induced res- 
onance). In this study we investigate the role of noise and 
diversity in a neuronal model of sleep-wake cycles based 
on the dynamics of the wake-promoting orexin neurons 
that is crucial for stability of wake and sleep states. We 
demonstrate that suitable levels of diversity introduced in 
the orexin neurons can significantly improve the quality 
of the sleep-wake cycle, and may be essential for proper 
sleep-wake periodicity. Noise, on the other hand, pro- 
vides only a mild improvement. 
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I. INTRODUCTION 

Disorder, which originates from both noise and diver- 
sity, is naturally present in all biological systems. In 
neuronal systems some examples are the random open- 
ing and closing of ion channels, the multitude of stochas- 
tic input currents in the neurons, and the diversity of 
shapes, sizes, and electrophysiological properties of the 
neurons [H, 0] ■ Disorder is often considered to be harmful 
to the systems' functioning and to information encoding. 
However, it was likewise repeatedly demonstrated that a 
certain level of disorder can facilitate signal encoding by 
enhancing system's response to an external stimuli. For 
instance, quenched diversity clearly shows its construc- 
tive role in the phenomenon of diversity-induced reso- 
nance, in which an assembly of heterogeneous excitable 
units presents an optimal response to an external forcing 
for a suitable intermediate degree of heterogeneity . 
Similar constructive effects can be observed in the pres- 
ence of noise. For example, interplay of noise and nonlin- 
ear forces produces the directed motion of motor proteins 
Q, order-disorder transitions, oscillations, and synchro- 
nization in assemblies of excitable units 043 > an d an op- 
timized system response in the ubiquitous phenomenon 
of stochastic resonance [13, [Hj], e.g. in ion-channels and 
neurons 

In the present study we examine the effects of noise and 
diversity (heterogeneity) in a physiologically based neu- 
ronal model of sleep-wake cycles fl8|. This model intro- 
duces a novel mechanism of the homeostatic regulation 
of sleep based on the dynamics of a wake-promoting neu- 
ropeptide orexin (also called hypocretin), assuming de- 
pression of orexinergic synapses during wakefulness and 
their recovery during sleep. This mechanism is based on 
the experimental findings of the essential role of orexin 
system in maintaining wakefulness and its ability to in- 
tegrate the sleep-wake relevant information coming from 
many brain areas [ljj [2(| and respond to changes in the 
body external and internal environments by encoding the 
body activity state, energy balance, sensory and emo- 
tional stimuli [HHll. 

In the original model interaction between only two 
representative neurons is simulated: the orexin neuron 
and the local glutamate neuron that are reciprocally con- 
nected to each other according to the experimentally es- 
tablished physiological connections [23j] . Both orexin and 
glutamate neurons are firing during wakefulness and are 
silent during sleep. The transitions between firing and si- 
lence are governed by the interplay between the circadian 
input and homeostatic mechanisms as initially proposed 
by Borbely [2^|. For simplicity, in this model only a sin- 
gle type of orexin neurotransmitter (instead of the two 
types actually known) is considered, and it is assumed 
that the system can be either in the wake state or in 
a generic non-Rapid Eye Movement sleep state, without 
specifying ultradian structure of sleep. Also this model 
did not consider noise effects, and diversity could not be 
included since there are only two neurons present. 



In the present paper we extend the above described 
two-neuron model to a more realistic multi-unit model 
with heterogeneous neurons. The aim of the study is to 
first of all investigate how the presence of diversity in the 
neuronal population affects sleep-wake transitions, since 
it is well-known that neurons are highly heterogeneous 
by their nature. In particular, within the orexin neurons 
population significant intrinsic diversity can be found: 
different electrophysiological properties, sizes in the di- 
ameter range 15-40 /im, and various shapes such as a 
spherical, fusiform, or multipolar [b| [25|] . Secondly, 
also stochastic fluctuations, representing current noise, 
are added to the model and the response of the system is 
studied for different levels of noise. The question natu- 
rally arises, to what extent noise and diversity are essen- 
tial ingredients for the functioning of assemblies of neu- 
rons and other complex systems, and what is the optimal 
level of noise and diversity required for the emergence of 
an optimal response to external stimuli. It is shown be- 
low that the model under study presents both diversity- 
induced resonance and stochastic resonance, but the for- 
mer appears more clear and robust, since it is always 
associated with a regular almost-periodic spiking-silence 
activity, rather than to the irregular random transitions 
characterizing the stochastic resonance regime. 



II. MATERIALS AND METHODS 

In this section the two-neuron model of sleep- wake cy- 
cles [TH is described and some examples of dynamics in 
the presence of an external periodic signal are illustrated. 
Further, this model is extended to account for multiple 
neurons dynamics and heterogeneity, and a simple quan- 
titative criterion to estimate the quality of a sleep-wake 
cycle is introduced. This criterion will be used in the 
Results section to compare sleep-wake cycles dynamics 
obtained at different parameter sets. 



A. The two-neuron model 

The original model of the homeostatic regulation of 
sleep has a minimal structure consisting of two represen- 
tative interacting neurons A and B, as depicted in Fig.[T] 
The neuron A simulates a representative neuron from 
the orexinergic neuronal population, while the neuron B 
represents a local glutamate interneuron (for details see 
jlH ) . The state of wakefulness or sleep is determined by 
the firing regime of neurons A and B, since these neurons 
are known to fire during wakefulness and be almost silent 
during sleep (see e.g. [22j|). 

Interaction between the neurons A and B takes place 
through glutamate and orexin neurotransmitters, as de- 
tailed below. The neuron A is acted upon by a stimu- 
lus in pace with the circadian rhythm, here treated as 
a periodic external signal — a simplification justified 
by its independence from the homeostatic process [26| . 
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FIG. 1: Scheme of the two-neuron model of the sleep- 
wake cycle |l8f] ■ The A — > B red arrow from the orexin- 
producing neuron A (red circle) to the neuron B (blue circle) 
represents the glutamate projection as well as the orexin pro- 
jection regulating the homeostatic process. The blue arrow 
represents the B — > A glutamate projection. The neuron A is 
also acted upon by a periodic signal representing the effect of 
the circadian clock. 



The homeostatic process itself is described by an addi- 
tional macroscopic variable M(t) simulating availability 
of orexin. 

Dynamics of the neurons A and B are based on a 
Hodgkin-Huxley-type model [27]. The membrane poten- 
tials of the neurons A (Va(£)) and B (Va(f)) are thus 
calculated as: 



C A 



dVA 
dt 

dV B 
3 dt 

9k[Vb 



= Ext + £a — Xa,L — ^A,Na — Xa,K — ^A,gl (1) 
= 4xt + CA - g^Wx-Ei] - S , Na[yA--E Na ]aA,Na 

- <?k [Va Ek]o,a,k ~ g s \[VA-E g \]a,A, s \, 

= Cb — Ib.L — ^B,Na — Ib,K — ^B,gl — Ex (2) 
= £b - 9l[Vb-E L } - 3Na[VB--E , Na]aB,Na 

Ek]o,b,k - 5 , g i[VB-£'gi]aB,gi - 5 , ox[VB--Box]aox, 



where C' p (p — A, B) are the membrane capacitances 
per unit area of the respective neurons, I a (a — 
L, Na, K, gl, ox) are the ionic currents, g a are the max- 
imum conductances, and E a are the equilibrium poten- 
tials. The capacitance values are taken as Ca — Cb = 
1/LtF/cm 2 . The values of all the other model parameters 
are listed in Table |TJ 

In the following we give a detailed explanation of dif- 
ferent parts of the model. 

• External forces. The current Ext acting on the neu- 
ron A and the noise currents £ p (t), p — A,B, can 
be considered as external forces, in the sense that 
they do not depend on the system variables. 

The external current Ext (t) is assumed to simulate 
a stimulus associated with the circadian rhythm. 
For simplicity in the present study a periodic pulse 
input is used to introduce circadian activation of 



the system: r, Ext(i) = Ext(i + T )- Such cur- 
rent can be interpreted as an awakening effect of 
an alarm clock or some other disturbance coming 
with a period of 24 hours. In the following we em- 
ploy a train of rectangular pulses with length To 
(to < t) and height I , as depicted in Fig.[2ftop, 



Iext(t) — Io j 
= 0, 



nr < t < jit + To , 

tit + T < t < (n + 1)t , (3) 



where n is an integer. This simple form is chosen 
because it is convenient for carrying out a system- 
atic study of the neuron response at different pa- 
rameters sets. However, it should be kept in mind 
that it represents a drastic simplification, and more 
realistic shapes of circadian currents can also be 
used [3. 

The noise term represents fluctuating cur- 

rents that are known to be always present in neu- 
rons. For simplicity, we assume zero-average Gaus- 
sian white-noise processes: 

(&>(*)> =0; 

(Ut)Zp'(s))=2D p S PiP ,S(t-s), p, P ' = A,B, (4) 

with Dp being the noise intensity. 

• Internal dynamics. The leakage, sodium, and 
potassium currents I p _ a (p = A,B; a — L, Na, K) 
in the equation of the neuron p depend only on the 
variables of the same neuron p and, thus, describe 
the neuronal internal dynamics. 

The leakage currents I p ^ = gh(V p — E^) repre- 
sent a flow of ions with a small conductance w 
0.1 fiS/cm 2 driving the membrane potential toward 
the negative value El ps — 60 mV. 

The depolarizing Na-currents E,,Na = <?Na(^p — 
J^Na) %>,Na have a maximum conductance g^n = 
3 //S/ cm 2 and a large positive equilibrium potential 
E-N& = 50 mV. The activation variables %,,Na(£)> 
with < a Pj Na < 1, represent the fraction of open 
ion-channels contributing to the Na current. Be- 
cause of their fast activation relative to the other 
time scales, the Na-current is assumed to be ac- 
tivated instantaneously, according to its voltage- 
dependency: 

Op,Na = $(S Na (V^ - WNa)) , (5) 

where <&(x) is the sigmoid function 

= 7T 1 ( \ > ( 6 ) 
1 + cxp(— X) 

is the steepness of the sigmoid function and 
W^Na is the half-activation potential. 

The repolarizing K-currents I Pt K — 9k(Yp — 
Ek) Qp.K are characterized by a maximum conduc- 
tance 3k — 4 /iS/cm 2 , a large negative equilib- 
rium potential Ek = — 90 mV, and a longer activa- 
tion time than the depolarizing Na-current, namely 



4 



TABLE I: Parameters of the two-neuron model [18]. 





Conductance 


Equilibrium 


Slope 


Threshold 


Time 






Potential 


Parameter 


Potential 


Scales 




(/iS/cm 2 ) 


(mV) 


(mV" 1 ) 


(mV) 


(ms) 


L (Leakage current) 


ffL = 0.1 


£ L = -60 








Na (Sodium current) 


<?Na = 3 


£ N a = 50 


S Na = 0.25 


W N a = -25 


(TNa « 0) 


K (Potassium current) 


<?K = 4 


E K = -90 


5 K = 0.25 


Wk = -25 


TK = 2 


gl (Glutamate current) 


ffgi = 0.15 


E g i = 50 


S g i = 1 


W g i = -20 


r g i = 30 


ox (Orexin current) 


gox = 0.135 
Sox = 0.2 


E ox = 50 


Sox. = 1 


W ox = -20 


r ox = 300 
r+ = 7500 
T~ = 920 


Periodic current 










r = 24000 
r ( , = 500 



tk = 2 ms. Consequently, the dynamics of the K- 
currents activation variables are modelled as 



da 



p,K 



(It 



1 



[ap,K ~ ®(S K {V P - W K ))] , (7) 



where $(x) is defined in Eq. ([6]). 

Couplings. The neurons A and B are mutually cou- 
pled by chemical synapses through the glutamate- 
induced (/ Plg i) and the orexin-induced (I ox ) cur- 
rents. Unlike the Na and K currents, J P)g i and 
I ox depend on the activity of both presynaptic 
and postsynaptic neurons. The activation variables 
a P! gi and a ox depend on the appearance of a spike 
in the presynaptic neuron, i.e. on the presynap- 
tic voltage. Additionally these currents depend on 
the voltage of the postsynaptic neuron, similarly to 
other ionic currents. Both glutamate and orexin 
are excitatory neurotransmitters, so they are as- 
sumed to open depolarizing ion channels, such as 
Na-channels. 

The activations of the glutamate-induced currents 
are modeled as: 

^ = --[a P ,zi-*(S sl (V p -W sl ))}, 
at T g i 

p = B if p= A; p = A if p = B. (8) 

This equation is similar to Eq. but has the 
important difference that the equilibrium value 
&(Sgi(Vp— W g i)) for the activation variable a P!g i de- 
pends on the membrane potential Vp of the other 
neuron p {p = B if p — A, p = A if p = B). The 
time constant r g i = 30 ms accounts for the delay 
coming from the activation of glutamate receptors, 
and the following activation of ion channels. 

The orexin-induced current represents the effect of 
orexin produced by the neuron A and acting on 
the neuron B. It is modeled in a form similar to the 
glutamate-induced current. This current provides 
a simplified description of the effects of orexin on 



the neuron B which appear after a complex series 
of processes, involving production of orexin in the 
soma of the neurons, its release in the synaptic cleft, 
and activation of G-protein coupled metabotropic 
receptors. The dynamics of the activation variable 
a ox (t) depend not only on the membrane poten- 
tial Va(£)j but are also related to the availability of 
orexin at time t, described by the additional vari- 
able M(t) (0 < M < 1). The dynamics of the 
variables a ox (i) and M(i) are defined by the equa- 
tions: 



da n 



1 



dt 
dM 
~~dl 



[Oox-M x $(Sox(VA-Wox))] , (9) 



_ _(M-1) 



—Mx<f>(S ox (V A -W ox )). 

7"ox 



(10) 



The term M x $(S ox (V A - W OK )) in the Eq. © 
reflects activation of the synaptic current due to 
appearance of a spike in the presynaptic neuron A. 
At the same time it determines the rate of orexin 
availability reduction in Eq. (TTU)) due to spiking of 
the neuron A with a time constant r^. The first 
term in Eq. (jlOl) determines recovery rate of the 
orexin availability with time constant r^. 

The meaning of the product M(t) x $(5 ox (Va.— 
W ox )) is that there is orexin-induced activity in 
neuron B if (1) there is enough orexin available 
above a critical threshold [M(t) > M cr i t i C ai], and 
(2) the neuron A is in the firing state [$(£) ~ 1]. 

The time constants accounting for the orexin 
dynamics are much longer than the time constants 
associated with ionic current terms. The time con- 
stant r ox of the homeostatic regulation process is 
even longer, being of the order of magnitude of the 
daily period r. 

For numerical convenience, simulations are made 
over rescaled daily and orexin time scales: the daily 
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FIG. 2: Response of the two-neuron model. Main variables and inter-spike times St p (p — A, B) versus time for a pulse 
height To = 0.895 mA (left) and Jo = 0.893 mA (right), see text for details. 



period was assumed to be r = 24 s, instead of 
t = 24 h, achieved through a suitable rescaling, 
which was applied to the orexin time scale r ox and 
the production and reduction times r^.. The other 
time parameters are left unchanged. Since such 



rescaled r ox and r^j are still much larger than any 
other time scale of the microscopic dynamics, the 
rescaling does not change the main results of the 
simulations. See [l8| for a detailed validation of 
such rescaling procedure. 
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All the parameter values for the currents are listed in 
Table HI It is assumed that the neurons A and B share the 
same parameter values, unless specified otherwise. Such 
an assumption is justified, because the major properties 
of these neurons required for the model are the tonic 
firing (periodic single spike activity) and silent states. 
Without any external input both neurons should be in a 
silent state, while they are brought to firing activity in 
response to depolarization. Therefore, change of param- 
eters in a physiologically allowed range would primarily 
lead to the different amount of depolarization needed to 
excite neurons, and would not affect the major outcomes 
of the simulations. 

The system defined above is essentially an excitable 
feedback system, i.e. both the external input of sufficient 
strength and the AB coupling are essential elements for 
maintaining firing activity of the neurons. Orexin-related 
dynamics, with the associated long time scales, are ex- 
pected to direct the homeostatic sleep process, which reg- 
ulates the sleep- wake transitions. The healthy sleep- wake 
cycles in this system are realized as follows: 

• Initiating wakefulness. A sufficiently strong or long 
external signal or a stimulus associated to the cir- 
cadian rhythm, e.g. the idealized rectangular pulse 
considered here, activates the system and induces 
firing activity in the neuron A. Due to the excita- 
tory synaptic connection from A to B, the neuron 
B is also activated. 

• Maintaining wakefulness. Once the pulse is finished 
and the external current is zero, the system remains 
in the wake state (i.e. both neurons A and B are 
firing) due to reciprocal excitation between the neu- 
rons. The firing activity lasts for a fraction q of the 
daily period r. Ideally one can assume q = 2/3, 
corresponding to 16 hours for a day of 24 hours, 
i.e. 16 seconds for the daily period r = 24 s with 
the time scales of the model considered here. 

• Initiating sleep. The firing stops in both neurons 
due to decreased availability of orexin according to 
the dynamics of M(t). This is associated witch the 
transition from wake to sleep. 

Two examples of the two-neuron model dynamics with- 
out noise are illustrated in Fig. [2] The left part of the 
figure represents the response obtained for a pulse length 
To = 500 ms and height Iq = 0.895 //A/cm 2 . In each 
period orexin is depleted during the neuronal activity 
and recovered while the neurons are silent. The stimu- 
lus parameters used in this example have been intention- 
ally chosen close to the critical firing threshold, so that 
by slightly reducing the pulse height or length, the peri- 
odic appearance of a continuous time interval of spiking 
regime is lost. Such case is demonstrated in the right 
hand side of the figure, where the current pulse height is 
slightly lower, Iq — 0.893 mA, while all the other param- 
eters are kept the same. There the prolonged wake state 



is induced only every other day, because the input is in- 
sufficient to induce sustained spiking at the same levels 
of orexin availability. By reducing the pulse amplitude 
or duration even further it is possible to observe different 
behaviors such as triple or higher-order periodicities. 




FIG. 3: Scheme of the heterogeneous model. Example 
of model system with Na = 5 orexin- producing neurons {Ai}, 
i — 1, . . . Na (red spheres) and one neuron B (blue sphere). 
The neurons A interact with each other through an all-to-all 
coupling (red lines). Blue and red projections have a meaning 
similar to those of Fig. [T] the neuron B is coupled to the 
neurons Ai through parallel glutamate projections, while each 
neuron Ai is coupled to neuron B through a glutamate and 
an orexin projection. The neurons A are also acted upon by 
a stimulus representing the effect of the circadian clock (gray 
arrows). 



B. The heterogeneous model 

As a step toward a more realistic model we generalize 
the two-neuron model into a heterogeneous multi-neuron 
model. For simplicity we first increase the number of 
orexin neurons only. To do this we replace the single 
neuron A by a set of Na neurons {A^} (i = 1, . . . , Na), 
while still maintaining only one neuron B. Also, in this 
paper we assume that the diversity is constant in time in 
order to consider the simplest case possible. 

In reality a certain level of heterogeneity is observed in 
all neuronal parameters. However, given that our model 
neurons are simple pacemaking neurons such diversifica- 
tion of different model parameters (in a physiologically 
allowed range) would simply lead to slightly different fir- 
ing rates of the neurons. This, in turn, will result in di- 
versity in activations of synaptic currents, which can be 
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mimicked by simply diversifying their activation thresh- 
olds. Thus, in the following we can limit ourselves to 
studying the effects of diversity in activation thresholds 
of synaptic currents without loss of generality. Further- 
more, as a first step, the heterogeneity is only introduced 
in the glutamate-induced currents to avoid having a too 
complicated system, which would become difficult to un- 
derstand. 

With regard to the coupling topology among the orexin 
neurons, so far there is no detailed experimental data. 
Therefore, for simplicity, we chose an all-to-all coupling 
via gap junctions, but other variations can be tested in 
the future. The intensity of the coupling has been chosen 
large enough to ensure that the neurons A, respond in 
pace to the external current. The equations of the two- 
neuron model are modified accordingly. 

• Dynamics of the neurons Aj. The membrane po- 
tentials V^\t) of the neurons Aj, are described by 
equations analogous to Eq. ([1]): 

cj<" 



(It 



= 1, 



<rW t\<-i i 
~ J A.L i J 



ext T~ q A 



(0 



A,L A,Na 



1 A,K J A : gl / 



-W TP 1„W .. .'I 



-9k[V^-E k ] a A z > - g g i[V A w -E gl ] 



l A,gl 



(0 



A 



v 



(11) 



The current terms are similar to those in the two- 
neuron model, apart from the additional coupling 
currents between two generic neurons Aj and A.,-, 
hi = k in t{V^-Vi j) ), with i,j = l,...,N A , where 
feint is the gap junctions conductance that can be 
treated as coupling strength. The currents' activa- 
tion variables a A N and ai'L are modeled in accord 
with the equations of the two-neuron model. Note 
that the specific values of the activation variables 
will be different for different neurons since they de- 
pend on voltages of each particular neuron Aj. 

For simplicity the same external current I cx t(t) 
given by Eq. ((3]) is assumed to act on all neurons 

Aj (see Fig. [3J . The noise terms £ A (*) as wen as 
the noise £b(£) acting on the neuron B (see be- 
low) are also defined similarly and assumed to be 
statistically independent from each other. For con- 
venience the properties of all stochastic forces are 
written together — 1, . . . , N A ): 

(d°(*)> = <&(*)> = o, 
(d i) (*)CB( S )> = o, 



(ei'(t)e A i, (s))=2D A s t , J s(t- S ), 

(&(t)tB(s)) = 2D B 5(t-s). 



(12) 



• Connections from the neuron B to the neurons Aj. 
The neuron B has glutamatergic synaptic inputs 
to each of the neurons Aj as depicted in Fig. [3] 
Diversity is introduced in the activation thresholds 
of the glutamate-induced currents according to the 
following equation for the activation variables: 



da 



0) 

A,gl 



dt 



1 



CO 

X A,gl 



The thresholds W B , adopt different values for each 
neuron Aj that are independently extracted from a 
probability distribution defined later in the text. 

• Connections from the neurons Aj to the neuron B. 
Each of the neurons Aj has synaptic projections to 
the neuron B. This is translated in the model by 
replacing the single glutamate- and orexin-induced 
currents with their averages such that Eqs. © and 
(El) for the activation variables become: 



Hs gi (v B -wi%)) 



(13) 



dav 



dt 

1 

r gi 
da B ,ox 



, N A 

«B. 8l -jE^#A (,) -<,)) 



,(14) 



dt 



1 



1 Na 

-^M«<D(S ox (y A W -iy ol )) 



N A 



j=i 



(15) 



Note that diversity is again introduced in the ac- 
tivation thresholds of the glutamate-induced cur- 
rents W^ l g i corresponding to heterogeneous (Aj— > 
B) synapses located at the neuron B. Due to the dif- 
ferences in the Aj neurons, the orexin availability 
function M'*'(t) is different for different neurons, 
although still following Eq. (fTU|) . 

The above described set of equations constitutes the 
multi-neuron heterogeneous model of the homeostatic 
regulation of sleep. Numerical results were obtained us- 
ing a variation of the Runge-Kutta 2nd-order method, 
which is suitable for equations with stochastic terms, 
namely the Heun method [28j . Identical initial condi- 
tions were assumed for all neurons, corresponding to a 
silent state. 



C. Quantifying the quality of the sleep-wake cycle 

In this section a heuristic criterion is introduced in or- 
der to evaluate and compare the quality of the system 
responses obtained for different external signals or inter- 
nal parameter values. 
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For this purpose, the period r is divided into a "day" 
wakefulness sub-period of length n = qr and a "night" 
sleep sub-period of length T2 = (1 — g)r, with r = n +T2. 
The quantity q is defined as a wafe fraction. A typical 
sleep- wake cycle with an eight-hour sleep sub-period has 
q = 2/3. For the day corresponding to the n-th period 
(nr,(n+l)r), the "day" is represented by the sub- interval 
(jit, nr + n) — (nr, (n + <7)t), which covers the first frac- 
tion g of the period, while the "night" extends in the 
complementary fraction (1 — q) of the period in the time 
interval (nT+T 2 ,(n + l)r) = ((n+g)r,(n+l)r). 

For each period n = 0, 1, . . . , we compute wakefulness 
time intervals At^ and Atn spent by the system in 
the wake state during the day, Atn , and night, Atn . 
The wake/sleep state is identified with the spiking/silent 
regime. A simple quantitative estimate of the quality of 
the sleep-wake cycle can, thus, be done through the fol- 
lowing linear function of the wakefulness time intervals, 



r(Ai«,A^) 



At« At^ 



n 



T-2 



(16) 



where At^ = J2n /N sp , a = 1,2, represent the 
average of the wakefulness time intervals during the day 
(q = 1) and during the night (a — 2), with N sp being 
the total number of periods of the simulation. 

The fractions At^ /r a (a = 1,2) can vary in the in- 
terval (0, 1); then the coefficient r in Eq. (fT6|) is limited 
in the interval (—1, 1). The maximum value r = 1 cor- 
responds to an optimal cycle with At^ = t± (wakeful- 
ness during the entire day) and At^ = (sleep during 
the entire night); any deviation from the optimal state 
(r = 1) comes either from values Ai^ 1 - 1 < t\ (implying 
some sleep during the day) or values At^ > (meaning 
at least some wakefulness in the night) . See the support- 
ing information in the Appendix for further details on 
the definition of the time intervals A^ 1 ), At^ and the 
coefficient r. 
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FIG. 4: Effect of noise in neurons {A^}. (A). Ten periods 
of the raster plots of neuron B for different intensities Da 
of the noise acting on neurons {Ai}. Vertical dashed lines 
mark the beginning of the pulses of the external current Text, 
see text for details. (B). Coefficient r, from Eq. ()16[). versus 
current noise intensity Da- 



III. RESULTS 



In this section we study how the presence of disor- 
der affects the system response and discuss the main dif- 
ferences compared to the two-neuron model. The term 
"disorder" is used here to refer to either noise, i.e. dis- 
order in time (stochastic terms in the external current), 
or diversity, i.e a quenched heterogeneity in the neuronal 
parameters. These two aspects are studied separately. 
For the sake of simplicity, we examine the response of 
the system to a periodic stimulus represented by a train 
of short rectangular pulses as defined in Materials and 
Methods. 

In each of the examples considered, the initial con- 
figuration in the absence of noise and diversity is the 
same as the sub-threshold state illustrated in Fig. [fright 
with a double-periodic response. It is obtained for a re- 
duced height of the current pulse Iq — 0.893 mA, while 



the other parameters are unchanged as given in Table 
|IJ The reason for starting from such an under-threshold 
non-optimal configuration is that it is most sensitive and, 
thus, best illustrates the effects of added noise or hetero- 
geneity. While a response with a double-periodicity may 
seem unrealistic, this starting configuration is intended 
to be an example of non-optimal response rather than a 
standard reference state. In fact, in realistic situations 
noise and heterogeneity are always present so that such 
a state without noise or diversity represents a hypothet- 
ical system that would be obtained if one could switch 
off noise or replace heterogeneous synapses with perfectly 
identical ones. The results presented below suggest that 
a multi-periodic sleep-wake cycle can be turned into a 
regular (single-periodic) one by adding a suitable degree 
of disorder. 
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FIG. 5: Effect of noise in neurons {A;}. Sample of four periods of some relevant variables and inter-spike times St p 
(p — A, B) versus time of neuron Ai and neuron B for an intensity of noise in neurons A Da = 1 mA (left) and Da = 5 mA 
(right). Compare Fig. [4] and see text for details. 



Effect of noise 



Noise in the neurons Ai 



Here we investigate the effects of the noise currents in 
the equations for the membrane potentials. For clarity 
only the cases in which noise currents are present either 
in the neurons A, or in the neuron B are considered. 



To study the effects of the noise currents ^(t), i = 
1, . . . , A^a acting only on the neurons A, (as per Eq. (ITTj) ) 
we set Db = 0. Also, no diversity in the characteristic 
parameters of neurons Aj is introduced. We have sim- 
ulated a system with Na — 20 identical neurons and a 
single B neuron on a time interval t € (0, 100 r). A raster 
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plot for the activity of the neuron B at different values 
of Da (indicated on the left) is shown in Fig. @JA. The 
plot shows that 

• for small Da ~ the system's configuration corre- 
sponds to the assumed non-optimal double-periodic 
solution; 

• the system's response becomes slightly more regu- 
lar and periodic as Da is increased, despite the fact 
that the neuron cannot initiate a firing event at the 
beginning of each period; 

• as Da becomes even larger the neuron B keeps fir- 
ing tonically for a longer and longer time interval 
(even longer than a single period) thus deteriorat- 
ing the general quality of the response. 

A sample of time dependence of the main variables in 
the interval (0,4r) for D A = 1 mV and Da = 2 mV 
is illustrated in Fig. [5j In general, the type of variabil- 
ity induced by noise currents acting on the neurons Aj 
affects both the firing initiation and, especially, its du- 
ration. However, it is difficult to establish an actual im- 
provement of the quality of such a response as a function 
of the noise intensity Da, as even the coefficient r, shown 
in Fig. 0}B, suggests only a mild stochastic resonant be- 
havior characterized by a wide plateau at intermediate 
values of Da- 



2. Noise in the neuron B 

Here we consider the complementary case, in which 
Da = and a current noise only affects the neuron B. A 
sample of raster plots of the membrane potential of the 
neuron B is depicted in Fig. [6jA in the time window t S 
(0, 10 r) for the values of the noise intensity Db indicated 
on the vertical axis. 

The raster plot in Fig. [BJ-A indicates that: 

• the smallest values of noise intensity Db ~ cor- 
respond to the double-periodic configuration dis- 
cussed above; 

• the response becomes periodic, and the length of 
the firing periods more regular for higher values of 

D B ; 

• at larger values of Db the state of sleep is frequently 
interrupted by almost isolated spikes at random 
times. 

A representative example of time dependence of selected 
variables of the neurons Ai and B are shown in Fig. [7] 
Note the different type of behavior induced by a high lev- 
els of noise acting on the neuron B, compared to the case 
in which noise acts on the neurons A^. In the first case 
irregular switching between the firing and silent states is 
observed more often, especially considering the transient 
firings in the otherwise silent sleep state. Furthermore, 
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FIG. 6: Effect of noise in neuron B. (A). Sample of ten 
periods of the raster plots of neuron B for different values of 
the intensity Db of the noise acting on neuron B. Vertical 
dashed lines mark the beginning of the pulses of the external 
current J e xt, see text for details. (B). Quality of the sleep- 
wake cycle from the coefficient r, Eq. (I16[) . versus current 
noise intensity Db- 



this random firing appears only in the neuron B, but is 
insufficient to also induce spiking in the Ai neuron. This 
activity may represents intermittent awakenings, which 
are likely due to the ability of noise to favor the igni- 
tion of spiking events. Such random spikes are not ob- 
served when noise acts on the neurons A^ only, even at 
much larger noise intensities. This may be related to the 
coupling between the neurons A^, which constrains them 
in the same (spiking or silent) state. In order to excite 
all neurons A^ together one would need an input signal 
affecting all of them in the same way, which is highly 
improbable in a realistic system. 

The dependence of the coefficient r on Db is shown 
in Fig. [BJ-B. Again, only a mild stochastic resonance be- 
havior is suggested by the data when varying the noise 
intensity. It should be noted that in this particular con- 
figuration with noise acting only on the neuron B, the 
response of the neuron B does not depend on the num- 
ber Na of homogeneous neurons {Aj}, due to the equiv- 
alence to the configuration of the two-neuron model, as 
we have checked numerically. Thus, the plots of neuron 
Ai in Fig. [7] are representative of all other neurons A^. 
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FIG. 7: Effect of noise in neuron B. Sample of four periods for some relevant variables and inter-spike times 8t p (p = A, B) 
versus time for neuron Ai and neuron B for an intensity of the noise acting on neuron B Db = 1 mA (left) and Db = 2 mA 
(right). Compare Fig. [6] and see text for details. 



In fact, the external current I ex t{t) as well as the cou- 
pling currents are the same for each neuron Aj, which 
produces the same response. According to the equations 
of the heterogeneous model, the effective current acting 
on the neuron B is the arithmetic average of the currents 
coming from the various neurons {Aj} and, therefore, co- 
incides with that of any single neuron Aj . We use here a 
homogeneous multi-neuron generalized model only for a 
better comparison and consistency with the rest of this 



study. 

B. Effects of heterogeneity 

The effects introduced by a heterogeneity in the neu- 
rons are dramatic compared to the effects of noise. The 
corresponding improvement of the system response for 
suitable intermediate amounts of diversity can be de- 
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FIG. 8: Effect of diversity in the B — > A» synapses. 

(A). Sample of ten periods of the raster plots of neuron B for 
different heterogeneity levels JH^B.gi in the B— >Ai glutamate 
synapse thresholds, see text for details. (B). Quality of the 
sleep- wake cycle from the coefficient r, Eq. (|16[l . for various 
threshold diversities 5Wb k i. 



tected very clearly. This is the main result of this paper 
and it is illustrated in this section. Noiseless neurons are 
assumed for easier estimation of the heterogeneity effects 
{D A = D B =0). 

As in the study of noise described above, we carry out 
the study of diversity starting from the same configura- 
tion with a non-optimal double-periodic response to the 
external periodic stimulus, corresponding to a zero diver- 
sity (homogeneous system). Heterogeneity is then intro- 
duced in the glutamate-induced currents, either in the 
thresholds wjf^ regulating the response of the A,— >B 

synapses at the neuron B or in the thresholds W B ( ^ of 
the B— >Ai synapses at the neurons A. This is done by 
randomly extracting values W from a probability den- 
sity f p (W) and assigning them to the threshold param- 
eters Wpg\ (p = A,B). The probability density used 
here has a bell-shape f p (W) = P((W - W p , g i)/SW p , g {), 



where P(x) cx l/cosh(a;) , the quantity W y 



(W) 



represents the average value, while SW p ^ g i measures the 
dispersion of the distribution f p (W) around the aver- 
age value and is related to the standard deviation a p by 

SWr ~ 



P;gl = 7Tf7 p /vl2. For further details see the support- 
ing information in the Appendix. The width 5W Plg \ is 
assumed in the following as the measure of neuronal di- 



versity. In order to carry out meaningful comparisons 
with the homogeneous (two-neuron) model, the average 
values are set equal to the corresponding parameters of 
the homogeneous two-neuron model, 



dWWf p {W) = W PigU p = A,B. (17) 



The other parameters are unchanged compared to the 
two-neuron model, see Table HI 



1. Diversity in the B— >Aj synapses (neurons Ai) 

(i) 

Diversifying the potential thresholds ^ implies het- 
erogeneous glutamate synapses located at the neurons 
A^ see Eq. (fT3|) and Fig. [3] That is, each neuron Ai 
responds in a different way to the stimulation from the 
neuron B. Notice that this is a truly heterogeneous sys- 
tem which cannot be reduced to an effective two-neuron 
model — as in the case of heterogeneous synapses at neu- 
ron B considered in the next section. We studied a sys- 
tem with Na = 20 neurons Aj with diversified threshold 
parameters W^'i, i = 1,...,Na- The system dynam- 
ics were examined for different sets of thresholds {Wg*^} 
extracted from distributions /b(W) with different widths 
<5Ws, g i but always the same average value W^B,gi = Ws.gi- 

The resulting raster plots of the activity of the neuron 
B are shown in Fig. [8} A, and a sample of time dependen- 
cies for the neurons Ai and B is shown in Fig. |9] The 
existence of an optimal degree of diversity, corresponding 
to a value <5WB, g i approximately between 1 and 1.5mV, 
can be clearly seen both from Fig. [8]-A and from the 
dependence of the coefficient r on the diversity degree 
<5W B , g i, in Fig.[5}B. 

The underlying mechanism leading the system from 
the double- to the single-periodic response as diversity 
is increased can be interpreted following the prototype 
mechanical model of diversity-induced resonance intro- 
duced in Ref. Q . In this model a set of interacting os- 
cillators moving in a bistable potential is subjected to an 
external periodic force, which pushes the system toward 
the left and the right barrier alternately. If the oscil- 
lators are identical, i.e. they have the same parameter 
values corresponding to an under-threshold regime, then 
the system of oscillators cannot perform jumps on the 
other site of the barrier under the action of the applied 
periodic force. However, when the parameters are diver- 
sified (keeping constant the corresponding average value) 
some oscillators respond more promptly to the force and 
jump to the other side of the barrier, gradually pulling 
the rest of the system. In the present case, each neuron 
Aj corresponds to a nonlinear oscillator of the example, 
while the parameter which is diversified is the activation 
(i) 

thresholds W-g gi of the glutamate-induced currents. 

To show that this is the actual mechanism in action, 
Fig. [TU] (left) illustrates the response of the heteroge- 
neous system by depicting the time dependence of the 
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FIG. 9: Effect of diversity in the B— >Ai synapses. Sample of four periods of some relevant variables and inter-spike times 
8t p (p = A, B) versus time for neuron Ai and neuron B for different levels of the threshold diversity 8Wb,& = 1 mV (left) and 
8Wb, s i = 5 mV (right). Compare Fig. [8] and see text for details. 



glutamate activation variables a£ g \(t) of the neurons Aj, 
i = 1,5, 10, 15, 20, with different values of the thresholds 
W'g'gj, at the beginning of a new period in the presence 
of the periodic current pulse. In Fig. [10] also the raster 
plots for all neurons in the same time interval are shown. 
One can notice that the activation variables oj^giW be- 
have differently from each other. Those associated to the 
lowest values of the activation threshold (indicated by 



small i values) respond stronger to the current pulse than 
those with the highest values of the threshold (largest 
values of i). The system is observed to reach the spik- 
ing regime faster than in the homogeneous case, which is 
shown in the right part of Fig. [TOl through the comparison 
between the glutamate activation variable of the homoge- 
neous system, a Agl (t), and the average activation vari- 
able (<XA,gl(i)) = ^YT 1 Si a Agl(*) °f tne heterogeneous 
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FIG. 10: Effect of diversity in the B— kA.* synapses. Comparison between the responses of the heterogeneous system (left 
column) and homogeneous system (right column) in the first part of the time period t/r £ (1, 2) during the action of the 500 ms 
long current pulse starting at ti/r — 1 and ending at t%jT ~ 1.021. (A) and (B) (top panels). External current pulse. (C) 
and (D) (central panels). Behavior of some representative glutamate activation variables of the heterogeneous system, a A g i(t) 
for i = 1,5, 10, 15, 20 (panel (C)), and the (common) time dependence aA, g i(i) of the homogeneous system activation variables 
(panel (D), black continuous curve); in the latter figure also the average value (aA,gl(i)) = N^ 1 Si a i giW °f ^he heterogeneous 
system (dashed grey curve) is shown for comparison. (E) and (F) (bottom panels). Raster plots of all the neurons of the 
system. See text for further details. 



system. Eventually, a Agl (t) — > and the homogeneous 
system goes back to the silent state, while the average ac- 
tivation variables of the heterogeneous system (and their 

average (&AgiW)) continue to oscillate around positive 
values, signaling the stability of the reached firing state. 



2. Diversity in the A — synapses (neuron B) 



In order to study the effects of added heterogeneity in 

the glutamate synapses located at the neuron B, one has 

(i) 

to diversify the potential threshold parameters W A g \, see 
Eq. (fT4"|) . For this particular case, it is possible to simplify 
the model into a two-neuron model with a single effective 
AB coupling. This is possible because heterogeneity only 
enters Eq. (fT4]) . while other model equations reduce to 
the same equations in the case of identical neurons Aj, 
so that all neurons Aj behave in the same way. Thus, the 



effective glutamate-induced current to the neuron B is 
1 Na 



dWf(W)$(S gl (V A -W)), (18) 



where Va is the common value of the membrane poten- 
tials of the neurons Aj. Here f(W) is the probability den- 
til 

sity of the corresponding thresholds W = W A1 , which 
is assumed to be the same bell-shaped probability den- 
sity f(W) = P((W - W)/5W) as discussed above, with 
W = W A gl and 8W = 6W A gl . 

We now consider two limiting cases of the effective cur- 
rent given by Eq. (fig)) . In the limit SW <C S^ 1 , when 

diversity is very small on the scale S^ 1 , it can be as- 
sumed that the following approximation holds, f(W) — 
P{{W ~W)/5W) w S(W-W) and the integral ^ can 
be reduced to the homogeneous result, 



I e s(V) ps / dW5(W~ W)$(S el (V A -W)) 

= $(S gl (V A -W)) . (19) 
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FIG. 11: Effect of diversity in the A;— >B synapses. (A). 
Sample of ten periods of raster plots of neuron B for different 
heterogeneity levels SWa^i in the Aj— >B glutamate synapse 
thresholds, see text for details. (B). Coefficient r, Eq. (|16l) . 
for various degrees of diversity SW A ,gl- 



In the complementary limit of high diversity level, SW 3> 
<S~j , the smooth function $(S g i(V A —W)) can be approx- 
imated with Heaviside step functions &(W—V A ), and the 
effective current becomes 

/ C ff 00 « [ dWf{w)e(w-v A ) 

dW P{{W - W)/5W) Q(W-V A ) 

= mV A -W)/SW). (20) 

This follows from the form of the chosen distribution, 
P(x) cx l/cosh(.T) 2 = — 4d$(x)/dx. Thus, in both these 
limiting cases the effective current can be written in the 
form Jeff (V) « $(A(V A -W 7 )), with A = S g i for SW < S^ 1 
and A = SW^ 1 for SW 3> S^ 1 . It is interesting that, as 
we have checked by numerical integration of Eq. (IT8l) , the 
same analytical form also holds for intermediate values 
of SW and S^ 1 , so to a very good approximation the ef- 
fective current can be written as I c g(V) w $(A(Va— W)), 
where the parameter A depends on the ratio <$W/ S^ 1 and 

varies monotonously between S g \ and as SW/S^ 1 

varies between and oo. 

The system's response at different levels of heterogene- 
ity, <$WA,gi, is presented in Fig. ITHA through the raster 
plots for the neuron B. A sample of time dependence is 



shown in Fig. [121 while Fig. fTTlB shows the dependence 
of the coefficient r on the diversity level. In Fig. [TTJ it 
can be seen that for small values of the diversity SWa, s \ 
the response of the neuron B presents the double peri- 
odicity of the reference configuration. Single periodicity 
is recovered for higher levels of diversity. At even higher 
values of 5Wa, s \, the coefficient r begins to decrease. The 
points of the curve corresponding to the highest values of 
r suggest an optimal degree of diversity SWa,%\ ~ 1 mV. 

The main difference compared to the case in which 
noise intensity is varied, is that the response of the sys- 
tem remains more regular also at the highest levels of 
diversity considered, i.e. without random spikes appear- 
ing during the silent state and with a typical cycle well 
shared between a day and a night sub-period. 



IV. DISCUSSION 

In the present work we have introduced a heteroge- 
neous multi-neuron version of the previously developed 
physiologically motivated model of the homeostatic reg- 
ulation of sleep. The multi-neuron model is composed 
of a population of conductance-based orexin-producing 
neurons and a single representative glutamatergic neu- 
ron. In this model the glutamatergic and orexinergic 
neurons are undergoing transitions between firing and 
silence depending on the external circadian input and in- 
ternal homeostatic mechanisms. These transitions corre- 
spond to the transitions between wake (firing) and sleep 
(silence), with the homeostatic mechanism being depen- 
dent on the availability of orexin. 

The specific aim of this study was to explore the ef- 
fects of noise and diversity in the regulation of sleep- 
wake cycles in such a model. It is clear that diversity 
and noise are integral parts of all biological systems, in- 
cluding the orexinergic neuronal population in the lateral 
hypothalamus. However, the role of disorder, and espe- 
cially diversity, is rarely considered in the physiologically 
based mathematical models of sleep-related systems |29r - 
Hfl . To our knowledge, diversity had so far been included 
only in one such model, i.e. the model of interacting cir- 
cadian oscillators , and here we present another exam- 
ple of the constructive role of diversity in regulation of 
sleep. 

We have demonstrated the existence of a diversity- 
induced resonance, leading to a clear and strong im- 
provement of the quality of the sleep- wake cycles, at 
a physiologically justified intermediate level of diver- 
sity of the orexin-producing neurons. However, only a 
mild improvement was found with varying noise inten- 
sity (stochastic resonance phenomena). 

We have considered the simplest system with only 20 
heterogeneous orexin neurons and one local glutamate 
neuron. Also we have used a very simple all-to-all net- 
work topology for the connections among orexinergic 
neurons. However, it can be expected that constructive 
effects of diversity will be found also in other model con- 
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FIG. 12: Effect of diversity in the A;— >B synapses. Sample of four periods of the time dependence of some relevant system 
variables and inter-spike times St p (p = A, B) of neuron Ai and neuron B for a threshold diversity 5W a, s i = 1 mV (left) and 
Wa, s \ = 5 mV (right). Compare Fig. [TTIand see text for details. 



figurations. In the future, more realistic modifications 
of the model with a larger population of glutamatergic 
neurons and more sophisticated inter-populations con- 
nections should be considered. Furthermore, in the fu- 
ture studies interplay between noise and diversity should 
likewise be investigated, since in nature both types of 
disorder are normally present. 

The validity of the result obtained within this model 
may be more general, since diversity-induced resonance is 



known to take place for suitable values of the parameters 
in general networks of interacting (non-linear) oscillators. 
A question then naturally arises: whether the phenomena 
encountered here could also characterize other systems 
where there is a coupling between two very different time 
scales or, in other words, if homeostatically regulated bi- 
ological systems may take advantage from a suitable level 
of heterogeneity of their components. 
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cycle. Therefore, on the difference of At^\ we have com- 
puted Af( 2 \ the wakefulness period during night, includ- 
ing also short wakefulness events (isolated spikes), since 
they break the sleep period, thus worsening the quality 
of sleep. Each isolated spike is assumed to contribute a 
conventional time interval, assumed as r max for simplic- 
ity. 

The difference between the two mentioned ways to 
evaluate wakefulness periods is relevant at high level of 
noise, which produces many isolated spikes. 



Appendix A: Supporting information: further 
details on the numerical simulation 



1. Evaluation of the wake and sleep times 

The quality of a sleep-wake cycle was estimates 
through the coefficient r = Ai^/ri — Av- 2 ' /t^, where 
At^ a ' represents the wakefulness time interval during the 
day (a = 1) or night (a = 2), while T\ and Ti represent 
the length of day and night, respectively. This coefficient 
can vary between a minimum value r = — 1, representing 
an exchange of wakefulness and sleep between day and 
night, and a maximum value r = 1, corresponding to the 
optimal situation with wakefulness during the whole day 
time interval At 1 - 1 ' = t\ and continuous sleep during the 
whole night interval, A^ 2 ' = 0. For a more realistic esti- 
mate of the quality of the sleep- wake cycle the two wake- 
fulness time intervals At^ and A^ 2 ) were computed in 
slightly different ways, as described below. 

Estimating wakefulness during the day: At^\ 
In order to estimate the wakefulness during the day, when 
an individual is supposed to be awake, we have consid- 
ered the quality of wakefulness by computing A6- 1 ' from 
states characterized by an "effective wakefulness" in a 
tonic firing regime. On the other hand, isolated spikes, 
which break a sleep period, did not contribute to the 
wakefulness interval and were neglected: such isolated 
spikes can at most represent a fragmented wakefulness 
state, similar to that of a narcoleptic individual. A tonic 
spiking state was recognized checking if the inter-spike 



time was smaller than a suitable threshold r n 



in this 



case such an inter-spike interval was added to the total 



wakefulness time interval At 



(i) 

wake ' 



On the other hand, if 



the corresponding inter-spike time was larger than the 



threshold 



then the two corresponding spikes were 



considered to be isolated and that time interval was ne- 
glected. Notice that in general inter-spike times are not 
universal and vary with external parameters; for this rea- 
son a value r max = 100 ms was chosen heuristically, con- 
sidering the typical working conditions of the system. 



Extraction of parameters with a given 
distribution 



The different values of the glutamate thresholds {Wi} 
of neurons {A^}, which make neurons A heterogeneous, 
were extracted using the probability distribution 



f(W) 



1 



2SWcosh z [(W - W)/SW] 



(Al) 



Here W = J dWW f(W) is the average value and the 
parameter SW measures the dispersion of the distribu- 
tion around W. It is proportional to the standard devia- 
tion cr, namely, it is related to the variance according to 
a 2 = ((W - W) 2 ) = tt 2 6W 2 /12. The function (fMJ can 
be integrated to obtain the lower cumulative distribution 
function F(x), 



F(W) = ~{1 



■taab[(W -W)/SW]} 
1 

1 +exp[-2(IF- W)/5W] ' 



(A2) 



Simulations were made by using various sets of A-neuron 
thresholds {Wi} obtained by rescaling the values of a 
same set of thresholds by SW. This can be done for 
example by inverting the same uniform distribution of 
values {Fi} and using distributions F(W) corresponding 
to the different desired values of SW. For all the param- 
eter sets {Wi} used the average value (W) was the same. 
For a few set of parameters we have checked that equiva- 
lent results are obtained by first extracting randomly the 
parameters from the distribution F(x) with the desired 
SW and then averaging the final results obtained. In the 
latter case the set of values {Fi} are extracted randomly 
in the interval F{ € (0, 1). 



Estimating wakefulness during the night: Ai^ 2 '. 
As for the night is concerned, it is the quality of sleep 
which contributes to the overall quality of the sleep- wake 
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